!=======================================================================
!
!BOP
!
! !MODULE: ice_step_mod
!
! !DESCRIPTION:
!
!  Contains CICE component driver routines common to all drivers.
!
! !REVISION HISTORY:
!  SVN:$Id: $
!
!  authors Elizabeth C. Hunke, LANL
!          Philip W. Jones, LANL
!          William H. Lipscomb, LANL
!
! 2008 ECH: created module by moving subroutines from drivers/cice4/
!
! !INTERFACE:
!
      module ice_step_mod
!
! !USES:
!
      use ice_atmo
      use ice_calendar
      use ice_communicate
      use ice_diagnostics
      use ice_domain
      use ice_dyn_evp
      use ice_fileunits
      use ice_flux
      use ice_grid
      use ice_history
      use ice_restart
      use ice_itd
      use ice_kinds_mod
      use ice_mechred
      use ice_ocean
      use ice_orbital
      use ice_shortwave
      use ice_state
      use ice_therm_itd
      use ice_therm_vertical
      use ice_timers
      use ice_transport_driver
      use ice_transport_remap
      use perf_mod, only: t_startf, t_stopf, t_barrierf

      implicit none
      private
      save

! !PUBLIC MEMBER FUNCTIONS:

      public :: step_therm2, step_dynamics, &
                prep_radiation, step_radiation
!
!EOP
!
!=======================================================================

      contains

!=======================================================================
!BOP
!
! !ROUTINE: prep_radiation - step pre-thermo radiation
!
! !DESCRIPTION: 
!
! !REVISION HISTORY:
!
! authors: Mariana Vertenstein, NCAR
!
! !INTERFACE:

      subroutine prep_radiation(dt)
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
      real (kind=dbl_kind), intent(in) :: &
         dt      ! time step
!
!EOP
!
      integer (kind=int_kind) :: &
         i,j,n,iblk    ! block index

      if (calc_Tsfc) then

         !$OMP PARALLEL DO PRIVATE(iblk)
         do iblk = 1, nblocks
            call prep_radiation_iblk(dt, iblk)
         end do
         !$OMP END PARALLEL DO

      else    ! .not. calc_Tsfc

         ! Initialize for safety
         do iblk = 1, nblocks
         do n = 1, ncat
         do j = 1, ny_block
         do i = 1, nx_block
            fswsfcn(i,j,n,iblk) = c0
            fswintn(i,j,n,iblk) = c0
            fswthrun(i,j,n,iblk) = c0
         enddo   ! i
         enddo   ! j
         enddo   ! ncat
            Iswabsn(:,:,:,iblk) = c0
            Sswabsn(:,:,:,iblk) = c0
         enddo   ! iblk

      endif    ! calc_Tsfc

      end subroutine prep_radiation

!=======================================================================
!BOP
!
! !ROUTINE: prep_radiation_iblk - step pre-thermo radiation
!
! !DESCRIPTION:
!
! !REVISION HISTORY:
!
! authors: David A. Bailey, NCAR
!
! !INTERFACE:

      subroutine prep_radiation_iblk (dt, iblk)
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
      real (kind=dbl_kind), intent(in) :: &
         dt      ! time step

      integer (kind=int_kind), intent(in) :: &
         iblk ! block index
!
!EOP
!
      integer (kind=int_kind) :: &
         i, j, ij    , & ! horizontal indices
         ilo,ihi,jlo,jhi, & ! beginning and end of physical domain
         n           , & ! thickness category index
         il1, il2    , & ! ice layer indices for eice
         sl1, sl2        ! snow layer indices for esno

      integer (kind=int_kind) :: &
         icells          ! number of cells with aicen > puny

      integer (kind=int_kind), dimension(nx_block*ny_block) :: &
         indxi, indxj    ! indirect indices for cells with aicen > puny

      ! snow variables for Delta-Eddington shortwave
      real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
         fsn             ! snow horizontal fraction
      real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr) :: &
         rhosnwn     , & ! snow density (kg/m3)
         rsnwn           ! snow grain radius (micro-meters)

      ! pond variables for Delta-Eddington shortwave
      real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
         fpn         , & ! pond fraction
         hpn             ! pond depth (m)

      real (kind=dbl_kind) :: netsw, netsw_old, ar

      type (block) :: &
         this_block      ! block information for current block

      logical (kind=log_kind) :: &
         l_stop          ! if true, abort the model

      integer (kind=int_kind) :: &
         istop, jstop    ! indices of grid cell where model aborts 

      l_stop = .false.

      fswfac(:,:,iblk) = c1

      this_block = get_block(blocks_ice(iblk),iblk)         
      ilo = this_block%ilo
      ihi = this_block%ihi
      jlo = this_block%jlo
      jhi = this_block%jhi

      !-----------------------------------------------------------------
      ! Compute netsw scaling factor (new netsw / old netsw)
      !-----------------------------------------------------------------

      do j = jlo, jhi
      do i = ilo, ihi
         if (aice(i,j,iblk) > c0 .and. scale_factor(i,j,iblk) > puny) then
            netsw = swvdr(i,j,iblk)*(c1 - alvdr_gbm(i,j,iblk)) &
                  + swvdf(i,j,iblk)*(c1 - alvdf_gbm(i,j,iblk)) &
                  + swidr(i,j,iblk)*(c1 - alidr_gbm(i,j,iblk)) &
                  + swidf(i,j,iblk)*(c1 - alidf_gbm(i,j,iblk))
            scale_factor(i,j,iblk) = netsw / scale_factor(i,j,iblk)
         else
            scale_factor(i,j,iblk) = c1
         endif
         fswfac(i,j,iblk) = scale_factor(i,j,iblk) ! for history
      enddo               ! i
      enddo               ! j

      do n = 1, ncat

      !-----------------------------------------------------------------
      ! Identify cells with nonzero ice area
      !-----------------------------------------------------------------

         icells = 0
         do j = jlo, jhi
         do i = ilo, ihi
            if (aicen(i,j,n,iblk) > puny) then
               icells = icells + 1
               indxi(icells) = i
               indxj(icells) = j
            endif
         enddo               ! i
         enddo               ! j

      !-----------------------------------------------------------------
      ! Scale absorbed solar radiation for change in net shortwave
      !-----------------------------------------------------------------

         il1 = ilyr1(n)
         il2 = ilyrn(n)
         sl1 = slyr1(n)
         sl2 = slyrn(n)

         do ij = 1, icells
            i = indxi(ij)
            j = indxj(ij)

            fswsfcn(i,j,n,iblk) = scale_factor(i,j,iblk)*fswsfcn (i,j,n,iblk)
            fswintn(i,j,n,iblk) = scale_factor(i,j,iblk)*fswintn (i,j,n,iblk)
            fswthrun(i,j,n,iblk)= scale_factor(i,j,iblk)*fswthrun(i,j,n,iblk)
            Sswabsn(i,j,sl1:sl2,iblk) = &
                    scale_factor(i,j,iblk)*Sswabsn(i,j,sl1:sl2,iblk)
            Iswabsn(i,j,il1:il2,iblk) = &
                    scale_factor(i,j,iblk)*Iswabsn(i,j,il1:il2,iblk)
         enddo
      enddo                  ! ncat

      end subroutine prep_radiation_iblk

!=======================================================================
!BOP
!
! !ROUTINE: step_therm2 - step post-coupler thermodynamics
!
! !DESCRIPTION:
!
!-----------------------------------------------------------------------
! Wrapper for driver for thermodynamic changes not needed for coupling:
! transport in thickness space, lateral growth and melting. Needed for 
! introducing OpenMP threading more simply.
!
! !REVISION HISTORY:
!
! author: Mariana Vertenstein, NCAR
!
! !INTERFACE:
!
      subroutine step_therm2 (dt)
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
      real (kind=dbl_kind), intent(in) :: &
         dt      ! time step
!
!EOP
!
      integer (kind=int_kind) :: &
         iblk,  &  ! block index
         i, j

      call ice_timer_start(timer_column)  ! column physics
      call ice_timer_start(timer_thermo)  ! thermodynamics
!      call ice_timer_start(timer_tmp)  ! temporary timer

      !$OMP PARALLEL DO PRIVATE(iblk)
      do iblk = 1, nblocks
         call step_therm2_iblk(dt, iblk)
      end do
      !$OMP END PARALLEL DO

      !-------------------------------------------------------------------
      ! Ghost cell updates for state variables.
      !-------------------------------------------------------------------

      call ice_timer_start(timer_bound)
      call bound_state (aicen, trcrn, &
                        vicen, vsnon, &
                        eicen, esnon)
      call ice_timer_stop(timer_bound)

      !$OMP PARALLEL DO PRIVATE(iblk,i,j)
      do iblk = 1, nblocks

      !-----------------------------------------------------------------
      ! Aggregate the updated state variables. 
      !----------------------------------------------------------------- 
 
         call aggregate (nx_block,          ny_block,             &
                         aicen(:,:,:,iblk), trcrn(:,:,:,:,iblk),  &
                         vicen(:,:,:,iblk), vsnon(:,:,  :,iblk),  &
                         eicen(:,:,:,iblk), esnon(:,:,  :,iblk),  &
                         aice (:,:,  iblk), trcr (:,:,:,  iblk),  &
                         vice (:,:,  iblk), vsno (:,:,    iblk),  &
                         eice (:,:,  iblk), esno (:,:,    iblk),  &
                         aice0(:,:,  iblk), tmask(:,:,    iblk),  &
                         ntrcr,             trcr_depend) 


      !-----------------------------------------------------------------
      ! Compute thermodynamic area and volume tendencies.
      !-----------------------------------------------------------------

         do j = 1, ny_block
         do i = 1, nx_block
            daidtt(i,j,iblk) = (aice(i,j,iblk) - daidtt(i,j,iblk)) / dt
            dvidtt(i,j,iblk) = (vice(i,j,iblk) - dvidtt(i,j,iblk)) / dt
         enddo
         enddo

      enddo                     ! iblk
      !$OMP END PARALLEL DO

!      call ice_timer_stop(timer_tmp)  ! temporary timer
      call ice_timer_stop(timer_thermo)  ! column physics
      call ice_timer_stop(timer_column)  ! column physics

      end subroutine step_therm2

!=======================================================================
!BOP
!
! !ROUTINE: step_therm2_iblk - step post-coupler thermodynamics
!
! !DESCRIPTION:
!
!-----------------------------------------------------------------------
! Driver for thermodynamic changes not needed for coupling:
! transport in thickness space, lateral growth and melting.
!
! NOTE: Ocean fluxes are initialized here.
!
! !REVISION HISTORY:
!
! author: William H. Lipscomb, LANL
!
! !INTERFACE:

      subroutine step_therm2_iblk (dt, iblk)
!
! !USES:
!
!
! !INPUT/OUTPUT PARAMETERS:
!
      real (kind=dbl_kind), intent(in) :: &
         dt      ! time step

      integer (kind=int_kind), intent(in) :: &
         iblk ! block index
!
!EOP
!
!lipscomb - delete hicen later?
!      real (kind=dbl_kind), &
!         dimension (nx_block,ny_block,ncat,max_blocks) :: &
!         hicen           ! ice thickness (m)

      integer (kind=int_kind) :: &
         ilo,ihi,jlo,jhi, & ! beginning and end of physical domain
         i, j, n

      integer (kind=int_kind) :: &
         icells          ! number of ice/ocean cells 

      integer (kind=int_kind), dimension(nx_block*ny_block) :: &
         indxi, indxj    ! indirect indices for ice/ocean cells

      type (block) :: &
         this_block      ! block information for current block

      logical (kind=log_kind) :: &
         l_stop          ! if true, abort model

      integer (kind=int_kind) :: &
         istop, jstop    ! indices of grid cell where model aborts

      l_stop = .false.

      this_block = get_block(blocks_ice(iblk),iblk)         
      ilo = this_block%ilo
      ihi = this_block%ihi
      jlo = this_block%jlo
      jhi = this_block%jhi

      !-----------------------------------------------------------------
      ! Let rain drain through to the ocean.
      !-----------------------------------------------------------------

      do j = 1, ny_block
      do i = 1, nx_block
         fresh     (i,j,iblk) = fresh(i,j,iblk)       &
              + frain(i,j,iblk)*aice(i,j,iblk)
      enddo
      enddo

      !-----------------------------------------------------------------
      ! Given thermodynamic growth rates, transport ice between
      ! thickness categories.
      !-----------------------------------------------------------------

      call ice_timer_start(timer_catconv, iblk) ! category conversions

      if (kitd == 1) then
      !-----------------------------------------------------------------
      ! Compute fractional ice area in each grid cell.
      !-----------------------------------------------------------------
         call aggregate_area (nx_block,          ny_block, &
                              aicen(:,:,:,iblk),           &
                              aice (:,:,  iblk), aice0(:,:,iblk))

      !-----------------------------------------------------------------
      ! Identify grid cells with ice.
      !-----------------------------------------------------------------

         icells = 0
         do j = jlo,jhi
         do i = ilo,ihi
            if (aice(i,j,iblk) > puny) then
               icells = icells + 1
               indxi(icells) = i
               indxj(icells) = j
            endif
         enddo
         enddo

         if (icells > 0) then

            call linear_itd (nx_block, ny_block,       &
                             icells, indxi, indxj,     &
                             ntrcr,  trcr_depend,      &
                             aicen_init(:,:,:,iblk),   &
                             vicen_init(:,:,:,iblk),   &
                             aicen     (:,:,:,iblk),   &
                             trcrn     (:,:,:,:,iblk), & 
                             vicen     (:,:,:,iblk),   &
                             vsnon     (:,:,:,iblk),   &
                             eicen     (:,:,:,iblk),   &
                             esnon     (:,:,:,iblk),   &
                             aice      (:,:,  iblk),   &
                             aice0     (:,:,  iblk),   &
                             l_stop,                   &
                             istop,    jstop)

            if (l_stop) then
               write (nu_diag,*) 'istep1, my_task, iblk =', &
                                  istep1, my_task, iblk
               write (nu_diag,*) 'Global block:', this_block%block_id
               if (istop > 0 .and. jstop > 0) &
                    write(nu_diag,*) 'Global i and j:', &
                                      this_block%i_glob(istop), &
                                      this_block%j_glob(jstop) 
               call abort_ice ('ice: Linear ITD error')
            endif

         endif

      endif  ! kitd

      call ice_timer_stop(timer_catconv, iblk)  ! category conversions

      !-----------------------------------------------------------------
      ! Add frazil ice growing in leads.
      !-----------------------------------------------------------------

      ! identify ice-ocean cells
      icells = 0
      do j = 1, ny_block
      do i = 1, nx_block
         if (tmask(i,j,iblk)) then
            icells = icells + 1
            indxi(icells) = i
            indxj(icells) = j
         endif
      enddo               ! i
      enddo               ! j

      call add_new_ice (nx_block,              ny_block, &
                        ntrcr,                 icells,   &
                        indxi,                 indxj,    &
                        tmask    (:,:,  iblk), dt,       &
                        aicen    (:,:,:,iblk),           &
                        trcrn    (:,:,:,:,iblk),         &
                        vicen    (:,:,:,iblk),           &
                        eicen    (:,:,:,iblk),           &
                        aice0    (:,:,  iblk),           &
                        aice     (:,:,  iblk),           &
                        frzmlt   (:,:,  iblk),           &
                        frazil   (:,:,  iblk),           &
                        frz_onset(:,:,  iblk), yday,     &
                        fresh    (:,:,  iblk),           &
                        fsalt    (:,:,  iblk),           &
                        Tf       (:,:,  iblk), l_stop,   &
                        istop, jstop)

      if (l_stop) then
         write (nu_diag,*) 'istep1, my_task, iblk =', &
                            istep1, my_task, iblk
         write (nu_diag,*) 'Global block:', this_block%block_id
         if (istop > 0 .and. jstop > 0) &
              write(nu_diag,*) 'Global i and j:', &
                                this_block%i_glob(istop), &
                                this_block%j_glob(jstop) 
         call abort_ice ('ice: add_new_ice error')
      endif

      !-----------------------------------------------------------------
      ! Melt ice laterally.
      !-----------------------------------------------------------------
      call lateral_melt (nx_block, ny_block,     &
                         ilo, ihi, jlo, jhi,     &
                         dt,                     &
                         fresh     (:,:,  iblk), &
                         fsalt     (:,:,  iblk), &    
                         fhocn     (:,:,  iblk), &
                         fsoot     (:,:,:,iblk), &
                         rside     (:,:,  iblk), &
                         meltl     (:,:,  iblk), &
                         aicen     (:,:,:,iblk), &
                         vicen     (:,:,:,iblk), &
                         vsnon     (:,:,:,iblk), &
                         eicen     (:,:,:,iblk), &
                         esnon     (:,:,:,iblk), &
                         trcrn     (:,:,:,:,iblk) )

      !-----------------------------------------------------------------
      ! For the special case of a single category, adjust the area and
      ! volume (assuming that half the volume change decreases the
      ! thickness, and the other half decreases the area).  
      !-----------------------------------------------------------------

!NOTE - this does not work - hicen_init is not defined - ECH

!         if (ncat==1) &
!              call reduce_area (nx_block, ny_block,     &
!                                ilo, ihi, jlo, jhi,     &
!                                tmask     (:,:,  iblk), &
!                                aicen     (:,:,:,iblk), &
!                                vicen     (:,:,:,iblk), &
!                                hicen_init(:,:,1,iblk), &
!                                hicen     (:,:,1,iblk)) 

      !-----------------------------------------------------------------
      ! ITD cleanup: Rebin thickness categories if necessary, and remove
      !  categories with very small areas.
      !-----------------------------------------------------------------

      call cleanup_itd (nx_block,             ny_block,             &
                        ilo, ihi,             jlo, jhi,             &
                        dt,                   ntrcr,                &
                        aicen   (:,:,:,iblk), trcrn (:,:,:,:,iblk), &
                        vicen   (:,:,:,iblk), vsnon (:,:,  :,iblk), &
                        eicen   (:,:,:,iblk), esnon (:,:,  :,iblk), &
                        aice0   (:,:,  iblk), aice      (:,:,iblk), &
                        trcr_depend,                                &
                        fresh   (:,:,  iblk), fsalt   (:,:,  iblk), &
                        fhocn   (:,:,  iblk), fsoot   (:,:,:,iblk), &
                        tr_aero,                                    &
                        heat_capacity,        l_stop,               &
                        istop,                jstop)

      if (l_stop) then
         write (nu_diag,*) 'istep1, my_task, iblk =', &
                            istep1, my_task, iblk
         write (nu_diag,*) 'Global block:', this_block%block_id
         if (istop > 0 .and. jstop > 0) &
              write(nu_diag,*) 'Global i and j:', &
                                this_block%i_glob(istop), &
                                this_block%j_glob(jstop) 
         call abort_ice ('ice: ITD cleanup error')
      endif

      end subroutine step_therm2_iblk

!=======================================================================
!BOP
!
! !ROUTINE: step_dynamics - step ice dynamics, transport, and ridging
!
! !DESCRIPTION:
!
! Run one time step of dynamics, horizontal transport, and ridging.
! NOTE: The evp and transport modules include boundary updates, so
!       they cannot be done inside a single block loop.  Ridging
!       and cleanup, on the other hand, are single-column operations. 
!       They are called with argument lists inside block loops
!       to increase modularity.
!
! !REVISION HISTORY:
!
! authors: William H. Lipscomb, LANL
!
! !INTERFACE:

      subroutine step_dynamics(dt_dyn,dt_thm)
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
      real (kind=dbl_kind), intent(in) :: &
         dt_dyn  , & ! dynamic time step
         dt_thm      ! thermodynamic time step for diagnostics
!
!EOP
!
      type (block) :: &
         this_block      ! block information for current block

      integer (kind=int_kind) :: & 
         iblk        , & ! block index 
         i,j         , & ! horizontal indices
         ilo,ihi,jlo,jhi ! beginning and end of physical domain

      integer (kind=int_kind) :: &
         icells          ! number of cells with aicen > puny

      integer (kind=int_kind), dimension(nx_block*ny_block) :: &
         indxi, indxj    ! indirect indices for cells with aicen > puny

      logical (kind=log_kind) :: &
         l_stop          ! if true, abort model

      integer (kind=int_kind) :: &
         istop, jstop    ! indices of grid cell where model aborts

      call init_history_dyn     ! initialize dynamic history variables

      dynCnt = dynCnt + 1

      !-----------------------------------------------------------------
      ! Elastic-viscous-plastic ice dynamics
      !-----------------------------------------------------------------

      call t_barrierf ('cice_dyn_evp_BARRIER',MPI_COMM_ICE)
      call t_startf ('cice_dyn_evp')
      if (kdyn == 1) call evp (dt_dyn)
      call t_stopf ('cice_dyn_evp')

      !-----------------------------------------------------------------
      ! Horizontal ice transport
      !-----------------------------------------------------------------

      call t_barrierf ('cice_dyn_horz_transport_BARRIER',MPI_COMM_ICE)
      call t_startf ('cice_dyn_horz_transport')
      if (advection == 'upwind') then
         call transport_upwind (dt_dyn)    ! upwind
      else
         call transport_remap (dt_dyn)     ! incremental remapping
      endif
      call t_stopf ('cice_dyn_horz_transport')

      !-----------------------------------------------------------------
      ! Ridging
      !-----------------------------------------------------------------

      call ice_timer_start(timer_column)
      call ice_timer_start(timer_ridge)
      call t_barrierf ('cice_dyn_ridge_BARRIER',MPI_COMM_ICE)
      call t_startf ('cice_dyn_ridge')

      l_stop = .false.

      !$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,&
      !$OMP	                icells,indxi,indxj,l_stop,istop,jstop)
      do iblk = 1, nblocks
         this_block = get_block(blocks_ice(iblk), iblk)
         ilo = this_block%ilo
         ihi = this_block%ihi
         jlo = this_block%jlo
         jhi = this_block%jhi

      !-----------------------------------------------------------------
      ! Identify ice-ocean cells.
      ! Note:  We can not define icells here using aice>puny because
      !        aice has not yet been updated since the transport (and
      !        it may be out of whack, which the ridging helps fix).-ECH
      !-----------------------------------------------------------------
           
         icells = 0
         do j = jlo, jhi
         do i = ilo, ihi
!PERF  There is a performance advantage here to only perform the ridging over a limited set of points
!PERF            if (tmask(i,j,iblk) .and. aice(i,j,iblk) > c0) then
            if (tmask(i,j,iblk)) then
               icells = icells + 1
               indxi(icells) = i
               indxj(icells) = j
            endif
         enddo               ! i
         enddo               ! j

         if (icells > 0) then

         call ridge_ice (nx_block,             ny_block,                 &
                         dt_dyn,               dt_thm,                   &
                         ntrcr,                icells,                   &
                         indxi,                indxj,                    &
!!                         Delt    (:,:,  iblk), divu      (:,:,  iblk), &
                         rdg_conv(:,:,  iblk), rdg_shear (:,:,  iblk),   &
                         aicen   (:,:,:,iblk), trcrn     (:,:,:,:,iblk), &
                         vicen   (:,:,:,iblk), vsnon     (:,:,:,iblk),   &
                         eicen   (:,:,:,iblk), esnon     (:,:,:,iblk),   &
                         aice0   (:,:,  iblk),                           &
                         trcr_depend,          l_stop,                   &
                         istop,                jstop,                    &   
                         dardg1dt(:,:,iblk),   dardg2dt  (:,:,iblk),     &
                         dvirdgdt(:,:,iblk),   opening   (:,:,iblk),     &
                         fresh   (:,:,iblk),   fhocn     (:,:,iblk),     &
                         fsoot   (:,:,:,iblk))

         if (l_stop) then
            write (nu_diag,*) 'istep1, my_task, iblk =', &
                               istep1, my_task, iblk
            write (nu_diag,*) 'Global block:', this_block%block_id
            if (istop > 0 .and. jstop > 0) &
                 write(nu_diag,*) 'Global i and j:', &
                                  this_block%i_glob(istop), &
                                  this_block%j_glob(jstop) 
            call abort_ice ('ice: Ridging error')
         endif

         endif

      enddo                     ! iblk
      !$OMP END PARALLEL DO

      call ice_timer_stop(timer_ridge)
      call t_stopf ('cice_dyn_ridge')

      call t_barrierf ('cice_dyn_column_BARRIER',MPI_COMM_ICE)
      call t_startf ('cice_dyn_column')

      !$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,&
      !$OMP	                icells,indxi,indxj,l_stop,istop,jstop)
      do iblk = 1, nblocks
         this_block = get_block(blocks_ice(iblk), iblk)
         ilo = this_block%ilo
         ihi = this_block%ihi
         jlo = this_block%jlo
         jhi = this_block%jhi

      !-----------------------------------------------------------------
      ! ITD cleanup: Rebin thickness categories if necessary, and remove
      !  categories with very small areas.
      !-----------------------------------------------------------------

         call cleanup_itd (nx_block,             ny_block,             &
                           ilo, ihi,             jlo, jhi,             &
                           dt_thm,               ntrcr,                &
                           aicen   (:,:,:,iblk), trcrn (:,:,:,:,iblk), &
                           vicen   (:,:,:,iblk), vsnon (:,:,  :,iblk), &
                           eicen   (:,:,:,iblk), esnon (:,:,  :,iblk), &
                           aice0   (:,:,  iblk), aice      (:,:,iblk), &
                           trcr_depend,                                &
                           fresh   (:,:,  iblk), fsalt   (:,:,  iblk), &
                           fhocn   (:,:,  iblk), fsoot   (:,:,:,iblk), &
                           tr_aero,                                    &
                           heat_capacity,        l_stop,               &
                           istop,                jstop)

         if (l_stop) then
            write (nu_diag,*) 'istep1, my_task, iblk =', &
                               istep1, my_task, iblk
            write (nu_diag,*) 'Global block:', this_block%block_id
            if (istop > 0 .and. jstop > 0) &
                 write(nu_diag,*) 'Global i and j:', &
                                  this_block%i_glob(istop), &
                                  this_block%j_glob(jstop) 
            call abort_ice ('ice: ITD cleanup error')
         endif

      enddo              ! iblk
      !$OMP END PARALLEL DO

      call t_stopf ('cice_dyn_column')

      !-------------------------------------------------------------------
      ! Ghost cell updates for state variables.
      !-------------------------------------------------------------------

      call t_barrierf ('cice_dyn_bound_BARRIER',MPI_COMM_ICE)
      call t_startf ('cice_dyn_bound')
      call ice_timer_start(timer_bound)
      call bound_state (aicen, trcrn, &
                        vicen, vsnon, &
                        eicen, esnon)
      call ice_timer_stop(timer_bound)
      call t_stopf ('cice_dyn_bound')

      call t_barrierf ('cice_dyn_agg_BARRIER',MPI_COMM_ICE)
      call t_startf ('cice_dyn_agg')

      !$OMP PARALLEL DO PRIVATE(iblk,this_block,ilo,ihi,jlo,jhi,i,j)
      do iblk = 1, nblocks

      !-----------------------------------------------------------------
      ! Aggregate the updated state variables. 
      !----------------------------------------------------------------- 
 
         call aggregate (nx_block,          ny_block,             &
                         aicen(:,:,:,iblk), trcrn(:,:,:,:,iblk),  &
                         vicen(:,:,:,iblk), vsnon(:,:,  :,iblk),  &
                         eicen(:,:,:,iblk), esnon(:,:,  :,iblk),  &
                         aice (:,:,  iblk), trcr (:,:,:,  iblk),  &
                         vice (:,:,  iblk), vsno (:,:,    iblk),  &
                         eice (:,:,  iblk), esno (:,:,    iblk),  &
                         aice0(:,:,  iblk), tmask(:,:,    iblk),  &
                         ntrcr,             trcr_depend) 

      !-----------------------------------------------------------------
      ! Compute dynamic area and volume tendencies.
      !-----------------------------------------------------------------

         this_block = get_block(blocks_ice(iblk),iblk)
         ilo = this_block%ilo
         ihi = this_block%ihi
         jlo = this_block%jlo
         jhi = this_block%jhi

         do j = jlo,jhi
         do i = ilo,ihi
            dvidtd(i,j,iblk) = (vice(i,j,iblk) - dvidtd(i,j,iblk)) /dt_dyn
            daidtd(i,j,iblk) = (aice(i,j,iblk) - daidtd(i,j,iblk)) /dt_dyn
         enddo
         enddo

      enddo              ! iblk
      !$OMP END PARALLEL DO

      call t_stopf ('cice_dyn_agg')
      call ice_timer_stop(timer_column)

      end subroutine step_dynamics

!=======================================================================
!BOP
!
! !ROUTINE: step_radiation - step pre-coupler radiation
!
! !DESCRIPTION:
!
! !REVISION HISTORY:
!
! authors: Mariana Vertenstein, NCAR
!
! !INTERFACE:

      subroutine step_radiation(dt)
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
      real (kind=dbl_kind), intent(in) :: &
         dt      ! time step
!
!EOP
!
      integer (kind=int_kind) :: &
         i,j,n,iblk    ! block index


      alvdr(:,:,:) = c0
      alvdf(:,:,:) = c0
      alidr(:,:,:) = c0
      alidf(:,:,:) = c0
      Sswabsn(:,:,:,:) = c0
#ifdef AEROFRC
      dalvdr_noaero(:,:,:) = c0
      dalvdf_noaero(:,:,:) = c0
      dalidr_noaero(:,:,:) = c0
      dalidf_noaero(:,:,:) = c0
      dSswabsn_noaero(:,:,:,:) = c0
#endif
#ifdef PONDFRC
      dalvdr_nopond(:,:,:) = c0
      dalvdf_nopond(:,:,:) = c0
      dalidr_nopond(:,:,:) = c0
      dalidf_nopond(:,:,:) = c0
      dSswabsn_nopond(:,:,:,:) = c0
#endif
!      do iblk = 1,nblocks
!         alvdr(:,:,iblk) = c0
!         alvdf(:,:,iblk) = c0
!         alidr(:,:,iblk) = c0
!         alidf(:,:,iblk) = c0
!         Sswabsn(:,:,:,iblk) = c0
!      enddo

      if (calc_Tsfc) then

         !$OMP PARALLEL DO PRIVATE(iblk)
         do iblk = 1, nblocks
            call step_radiation_iblk(dt, iblk)
         end do
         !$OMP END PARALLEL DO

      else    ! .not. calc_Tsfc

         ! Initialize for safety
         do iblk = 1, nblocks
         do n = 1, ncat
         do j = 1, ny_block
         do i = 1, nx_block
            alvdrn(i,j,n,iblk) = c0
            alidrn(i,j,n,iblk) = c0
            alvdfn(i,j,n,iblk) = c0
            alidfn(i,j,n,iblk) = c0
            fswsfcn(i,j,n,iblk) = c0
            fswintn(i,j,n,iblk) = c0
            fswthrun(i,j,n,iblk) = c0
         enddo   ! i
         enddo   ! j
         enddo   ! ncat
            Iswabsn(:,:,:,iblk) = c0
            Sswabsn(:,:,:,iblk) = c0
         enddo   ! iblk

      endif    ! calc_Tsfc

      end subroutine step_radiation

!=======================================================================
!BOP
!
! !ROUTINE: step_radiation_iblk - step pre-coupler radiation
!
! !DESCRIPTION:
!
! !REVISION HISTORY:
!
! authors: David A. Bailey, NCAR
!
! !INTERFACE:

      subroutine step_radiation_iblk (dt, iblk)
!
! !USES:
!
! !INPUT/OUTPUT PARAMETERS:
!
      real (kind=dbl_kind), intent(in) :: &
         dt      ! time step

      integer (kind=int_kind), intent(in) :: &
         iblk ! block index
!
!EOP
!
      integer (kind=int_kind) :: &
         i, j, ij    , & ! horizontal indices
         ilo,ihi,jlo,jhi, & ! beginning and end of physical domain
         n           , & ! thickness category index
         il1, il2    , & ! ice layer indices for eice
         sl1, sl2        ! snow layer indices for esno

      integer (kind=int_kind) :: &
         icells          ! number of cells with aicen > puny

      integer (kind=int_kind), dimension(nx_block*ny_block) :: &
         indxi, indxj    ! indirect indices for cells with aicen > puny

      ! snow variables for Delta-Eddington shortwave
      real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
         fsn             ! snow horizontal fraction
      real (kind=dbl_kind), dimension (nx_block,ny_block,nslyr) :: &
         rhosnwn     , & ! snow density (kg/m3)
         rsnwn           ! snow grain radius (micro-meters)

      ! pond variables for Delta-Eddington shortwave
      real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
         fpn         , & ! pond fraction
         hpn             ! pond depth (m)

      type (block) :: &
         this_block      ! block information for current block

      logical (kind=log_kind) :: &
         l_stop          ! if true, abort the model

      integer (kind=int_kind) :: &
         istop, jstop    ! indices of grid cell where model aborts 

      l_stop = .false.

      this_block = get_block(blocks_ice(iblk),iblk)         
      ilo = this_block%ilo
      ihi = this_block%ihi
      jlo = this_block%jlo
      jhi = this_block%jhi

      !-----------------------------------------------------------------
      ! Compute cosine of solar zenith angle.
      ! This is used by the delta-Eddington shortwave module.
      ! Albedos are aggregated in merge_fluxes only for cells w/ coszen > 0.
      ! For basic shortwave, simply set coszen to a constant between 0 and 1.
      !-----------------------------------------------------------------

      if (trim(shortwave) == 'dEdd') then ! delta Eddington

         ! identify ice-ocean cells
         icells = 0
         do j = 1, ny_block
         do i = 1, nx_block
            if (tmask(i,j,iblk)) then
               icells = icells + 1
               indxi(icells) = i
               indxj(icells) = j
            endif
         enddo               ! i
         enddo               ! j

         call compute_coszen (nx_block,         ny_block,       &
                              icells,                           &
                              indxi,            indxj,          &
                              tlat  (:,:,iblk), tlon(:,:,iblk), &
                              coszen(:,:,iblk), dt)

      else                     ! basic (ccsm3) shortwave
         coszen(:,:,iblk) = p5 ! sun above the horizon
      endif

      do n = 1, ncat

      !-----------------------------------------------------------------
      ! Identify cells with nonzero ice area
      !-----------------------------------------------------------------
           
         icells = 0
         do j = jlo, jhi
         do i = ilo, ihi
            if (aicen(i,j,n,iblk) > puny) then
               icells = icells + 1
               indxi(icells) = i
               indxj(icells) = j
            endif
         enddo               ! i
         enddo               ! j

      !-----------------------------------------------------------------
      ! Solar radiation: albedo and absorbed shortwave
      !-----------------------------------------------------------------

         il1 = ilyr1(n)
         il2 = ilyrn(n)
         sl1 = slyr1(n)
         sl2 = slyrn(n)

         if (trim(shortwave) == 'dEdd') then   ! delta Eddington

      ! note that rhoswn, rsnw, fp, hp and Sswabs ARE NOT dimensioned with ncat
      ! BPB 19 Dec 2006

            ! set snow properties
            call shortwave_dEdd_set_snow(nx_block, ny_block,           &
                              icells,                                  &
                              indxi,               indxj,              &
                              aicen(:,:,n,iblk),   vsnon(:,:,n,iblk),  &
                              trcrn(:,:,nt_Tsfc,n,iblk), fsn,          &
                              rhosnwn,             rsnwn)


            if (.not. tr_pond) then

               ! set pond properties
               call shortwave_dEdd_set_pond(nx_block, ny_block,            &
                                 icells,                                   &
                                 indxi,               indxj,               &
                                 aicen(:,:,n,iblk),                        &
                                 trcrn(:,:,nt_Tsfc,n,iblk),                &
                                 fsn,                 fpn,                 &
                                 hpn)

            else

               fpn(:,:) = apondn(:,:,n,iblk)
               hpn(:,:) = hpondn(:,:,n,iblk)

            endif

#ifdef AEROFRC
            ! Diagnose aerosol forcing

            if (tr_aero) then

            tr_aero = .false.

            call shortwave_dEdd(nx_block,        ny_block,            &
                                icells,                                 &
                                indxi,             indxj,               &
                                coszen(:,:, iblk),                      &
                                aicen(:,:,n,iblk), vicen(:,:,n,iblk),   &
                                vsnon(:,:,n,iblk), fsn,                 &
                                rhosnwn,           rsnwn,               &
                                fpn,               hpn,                 &
                                trcrn(:,:,:,n,iblk), tarea(:,:,iblk),   &
                                swvdr(:,:,  iblk), swvdf(:,:,  iblk),   &
                                swidr(:,:,  iblk), swidf(:,:,  iblk),   &
                                dalvdrn_noaero(:,:,n,iblk),              &
                                dalvdfn_noaero(:,:,n,iblk),              &
                                dalidrn_noaero(:,:,n,iblk),              &
                                dalidfn_noaero(:,:,n,iblk),              &
                                dfswsfcn_noaero(:,:,n,iblk),             &
                                dfswintn_noaero(:,:,n,iblk),             &
                                dfswthrun_noaero(:,:,n,iblk),            &
                                dSswabsn_noaero(:,:,sl1:sl2,iblk),       &
                                dIswabsn_noaero(:,:,il1:il2,iblk),       &
                                dalbicen_noaero(:,:,n,iblk),             &
                                dalbsnon_noaero(:,:,n,iblk),             &
                                dalbpndn_noaero(:,:,n,iblk))

            tr_aero = .true.

            endif
#endif

#ifdef PONDFRC
            ! Diagnose pond forcing

            if (tr_pond) then

            tr_aero = .false.

            fpn(:,:) = c0
            hpn(:,:) = c0

            call shortwave_dEdd(nx_block,        ny_block,            &
                                icells,                                 &
                                indxi,             indxj,               &
                                coszen(:,:, iblk),                      &
                                aicen(:,:,n,iblk), vicen(:,:,n,iblk),   &
                                vsnon(:,:,n,iblk), fsn,                 &
                                rhosnwn,           rsnwn,               &
                                fpn,               hpn,                 &
                                trcrn(:,:,:,n,iblk), tarea(:,:,iblk),   &
                                swvdr(:,:,  iblk), swvdf(:,:,  iblk),   &
                                swidr(:,:,  iblk), swidf(:,:,  iblk),   &
                                dalvdrn_nopond(:,:,n,iblk),              &
                                dalvdfn_nopond(:,:,n,iblk),              &
                                dalidrn_nopond(:,:,n,iblk),              &
                                dalidfn_nopond(:,:,n,iblk),              &
                                dfswsfcn_nopond(:,:,n,iblk),             &
                                dfswintn_nopond(:,:,n,iblk),             &
                                dfswthrun_nopond(:,:,n,iblk),            &
                                dSswabsn_nopond(:,:,sl1:sl2,iblk),       &
                                dIswabsn_nopond(:,:,il1:il2,iblk),       &
                                dalbicen_nopond(:,:,n,iblk),             &
                                dalbsnon_nopond(:,:,n,iblk),             &
                                dalbpndn_nopond(:,:,n,iblk))

            fpn(:,:) = apondn(:,:,n,iblk)
            hpn(:,:) = hpondn(:,:,n,iblk)

            tr_aero = .true.

            endif
#endif
#ifdef CCSM3FRC
            call shortwave_ccsm3(nx_block,     ny_block,           &
                           icells,                                 &
                           indxi,             indxj,               &
                           aicen(:,:,n,iblk), vicen(:,:,n,iblk),   &
                           vsnon(:,:,n,iblk),                      &
                           trcrn(:,:,nt_Tsfc,n,iblk),              &
                           swvdr(:,:,  iblk), swvdf(:,:,  iblk),   &
                           swidr(:,:,  iblk), swidf(:,:,  iblk),   &
                           dalvdrn_ccsm3(:,:,n,iblk),               &
                           dalidrn_ccsm3(:,:,n,iblk),               &
                           dalvdfn_ccsm3(:,:,n,iblk),               &
                           dalidfn_ccsm3(:,:,n,iblk),               &
                           dfswsfcn_ccsm3(:,:,n,iblk),              &
                           dfswintn_ccsm3(:,:,n,iblk),              &
                           dfswthrun_ccsm3(:,:,n,iblk),             &
                           dIswabsn_ccsm3(:,:,il1:il2,iblk),        &
                           dalbicen_ccsm3(:,:,n,iblk),              &
                           dalbsnon_ccsm3(:,:,n,iblk))

#endif
            call shortwave_dEdd(nx_block,        ny_block,            &
                                icells,                                 &
                                indxi,             indxj,               &
                                coszen(:,:, iblk),                      &
                                aicen(:,:,n,iblk), vicen(:,:,n,iblk),   &
                                vsnon(:,:,n,iblk), fsn,                 &
                                rhosnwn,           rsnwn,               &
                                fpn,               hpn,                 &
                                trcrn(:,:,:,n,iblk), tarea(:,:,iblk),   &
                                swvdr(:,:,  iblk), swvdf(:,:,  iblk),   &
                                swidr(:,:,  iblk), swidf(:,:,  iblk),   &
                                alvdrn(:,:,n,iblk),alvdfn(:,:,n,iblk),  &
                                alidrn(:,:,n,iblk),alidfn(:,:,n,iblk),  &
                                fswsfcn(:,:,n,iblk),fswintn(:,:,n,iblk),&
                                fswthrun(:,:,n,iblk),                   &
                                Sswabsn(:,:,sl1:sl2,iblk),              &
                                Iswabsn(:,:,il1:il2,iblk),              &
                                albicen(:,:,n,iblk),                    &
                                albsnon(:,:,n,iblk),albpndn(:,:,n,iblk))

         else

            Sswabsn(:,:,sl1:sl2,iblk) = c0

            call shortwave_ccsm3(nx_block,     ny_block,           &
                           icells,                                 &
                           indxi,             indxj,               &
                           aicen(:,:,n,iblk), vicen(:,:,n,iblk),   &
                           vsnon(:,:,n,iblk),                      &
                           trcrn(:,:,nt_Tsfc,n,iblk),              &
                           swvdr(:,:,  iblk), swvdf(:,:,  iblk),   &
                           swidr(:,:,  iblk), swidf(:,:,  iblk),   &
                           alvdrn(:,:,n,iblk),alidrn(:,:,n,iblk),  &
                           alvdfn(:,:,n,iblk),alidfn(:,:,n,iblk),  &
                           fswsfcn(:,:,n,iblk),fswintn(:,:,n,iblk),&
                           fswthrun(:,:,n,iblk),                   &
                           Iswabsn(:,:,il1:il2,iblk),              &
                           albicen(:,:,n,iblk),albsnon(:,:,n,iblk))

         endif

#ifdef AEROFRC
            dalvdrn_noaero(:,:,n,iblk) = dalvdrn_noaero(:,:,n,iblk)-alvdrn(:,:,n,iblk)
            dalvdfn_noaero(:,:,n,iblk) = dalvdfn_noaero(:,:,n,iblk)-alvdfn(:,:,n,iblk)
            dalidrn_noaero(:,:,n,iblk) = dalidrn_noaero(:,:,n,iblk)-alidrn(:,:,n,iblk)
            dalidfn_noaero(:,:,n,iblk) = dalidfn_noaero(:,:,n,iblk)-alidfn(:,:,n,iblk)
            dfswsfcn_noaero(:,:,n,iblk) = dfswsfcn_noaero(:,:,n,iblk)-fswsfcn(:,:,n,iblk)
            dfswintn_noaero(:,:,n,iblk) = dfswintn_noaero(:,:,n,iblk)-fswintn(:,:,n,iblk)
            dfswthrun_noaero(:,:,n,iblk) = dfswthrun_noaero(:,:,n,iblk)-fswthrun(:,:,n,iblk)
            dfswabsn_noaero(:,:,n,iblk) = dfswsfcn_noaero(:,:,n,iblk)+dfswintn_noaero(:,:,n,iblk)+dfswthrun_noaero(:,:,n,iblk)
            dalbicen_noaero(:,:,n,iblk) = dalbicen_noaero(:,:,n,iblk)-albicen(:,:,n,iblk)
            dalbsnon_noaero(:,:,n,iblk) = dalbsnon_noaero(:,:,n,iblk)-albsnon(:,:,n,iblk)
            dalbpndn_noaero(:,:,n,iblk) = dalbpndn_noaero(:,:,n,iblk)-albpndn(:,:,n,iblk)
            dSswabsn_noaero(:,:,sl1:sl2,iblk) = dSswabsn_noaero(:,:,sl1:sl2,iblk)-Sswabsn(:,:,sl1:sl2,iblk)
            dIswabsn_noaero(:,:,il1:il2,iblk) = dIswabsn_noaero(:,:,il1:il2,iblk)-Iswabsn(:,:,il1:il2,iblk)
#endif
#ifdef CCSM3FRC
            dalvdrn_ccsm3(:,:,n,iblk) = dalvdrn_ccsm3(:,:,n,iblk)-alvdrn(:,:,n,iblk)
            dalvdfn_ccsm3(:,:,n,iblk) = dalvdfn_ccsm3(:,:,n,iblk)-alvdfn(:,:,n,iblk)
            dalidrn_ccsm3(:,:,n,iblk) = dalidrn_ccsm3(:,:,n,iblk)-alidrn(:,:,n,iblk)
            dalidfn_ccsm3(:,:,n,iblk) = dalidfn_ccsm3(:,:,n,iblk)-alidfn(:,:,n,iblk)
            dfswsfcn_ccsm3(:,:,n,iblk) = dfswsfcn_ccsm3(:,:,n,iblk)-fswsfcn(:,:,n,iblk)
            dfswintn_ccsm3(:,:,n,iblk) = dfswintn_ccsm3(:,:,n,iblk)-fswintn(:,:,n,iblk)
            dfswthrun_ccsm3(:,:,n,iblk) = dfswthrun_ccsm3(:,:,n,iblk)-fswthrun(:,:,n,iblk)
            dfswabsn_ccsm3(:,:,n,iblk) = dfswsfcn_ccsm3(:,:,n,iblk)+dfswintn_ccsm3(:,:,n,iblk)+dfswthrun_ccsm3(:,:,n,iblk)
            dalbicen_ccsm3(:,:,n,iblk) = dalbicen_ccsm3(:,:,n,iblk)-albicen(:,:,n,iblk)
            dalbsnon_ccsm3(:,:,n,iblk) = dalbsnon_ccsm3(:,:,n,iblk)-albsnon(:,:,n,iblk)
            dIswabsn_ccsm3(:,:,il1:il2,iblk) = dIswabsn_ccsm3(:,:,il1:il2,iblk)-Iswabsn(:,:,il1:il2,iblk)
#endif
#ifdef PONDFRC
            dalvdrn_nopond(:,:,n,iblk) = dalvdrn_nopond(:,:,n,iblk)-alvdrn(:,:,n,iblk)
            dalvdfn_nopond(:,:,n,iblk) = dalvdfn_nopond(:,:,n,iblk)-alvdfn(:,:,n,iblk)
            dalidrn_nopond(:,:,n,iblk) = dalidrn_nopond(:,:,n,iblk)-alidrn(:,:,n,iblk)
            dalidfn_nopond(:,:,n,iblk) = dalidfn_nopond(:,:,n,iblk)-alidfn(:,:,n,iblk)
            dfswsfcn_nopond(:,:,n,iblk) = dfswsfcn_nopond(:,:,n,iblk)-fswsfcn(:,:,n,iblk)
            dfswintn_nopond(:,:,n,iblk) = dfswintn_nopond(:,:,n,iblk)-fswintn(:,:,n,iblk)
            dfswthrun_nopond(:,:,n,iblk) = dfswthrun_nopond(:,:,n,iblk)-fswthrun(:,:,n,iblk)
            dfswabsn_nopond(:,:,n,iblk) = dfswsfcn_nopond(:,:,n,iblk)+dfswintn_nopond(:,:,n,iblk)+dfswthrun_nopond(:,:,n,iblk)
            dalbicen_nopond(:,:,n,iblk) = dalbicen_nopond(:,:,n,iblk)-albicen(:,:,n,iblk)
            dalbsnon_nopond(:,:,n,iblk) = dalbsnon_nopond(:,:,n,iblk)-albsnon(:,:,n,iblk)
            dalbpndn_nopond(:,:,n,iblk) = dalbpndn_nopond(:,:,n,iblk)-albpndn(:,:,n,iblk)
            dSswabsn_nopond(:,:,sl1:sl2,iblk) = dSswabsn_nopond(:,:,sl1:sl2,iblk)-Sswabsn(:,:,sl1:sl2,iblk)
            dIswabsn_nopond(:,:,il1:il2,iblk) = dIswabsn_nopond(:,:,il1:il2,iblk)-Iswabsn(:,:,il1:il2,iblk)
#endif
      enddo                  ! ncat

      end subroutine step_radiation_iblk

!=======================================================================

      end module ice_step_mod

!=======================================================================
